#pragma once

#include "Mesh.h"
#include "Heap.h"
#include "Matrix.h"

/**
 * @file MeshSimplify.h
 * @author xingyifan
 * @date 2022-12-17 23:00
 *
 * @description: 网格简化算法
 */

using namespace Mesh;
using namespace LinearAlgebra;

// 获得顶点的二次误差矩阵
Matrix QEMatrix(vertex3f *v)
{
    Matrix Q(4, 4);
    // 获得顶点的相邻面，计算每个面的法向
    std::deque<face *> fList = adjacentFace(v);
    for (int j = 0; j < fList.size(); ++j)
    {
        // 要计算平面方程 ax+by+cz+d=0
        vector3f normal = faceNormal(fList[j]);
        double d = -dot(normal, {v->pos.x, v->pos.y, v->pos.z});
        Vector p = {normal.x, normal.y, normal.z, d};
        Q = Q + p * p.transpose();
    }
    return Q;
}

// 曲面简化算法 QEM，不考虑非流形网格：有效对总是边
bool simplify(HEMesh *mesh, int times)
{
    // 首先要预处理获得所有有效对。由于不考虑非流形，有效对就是所有半边，什么也不用做
    // 然后对每个顶点计算一个 4x4 误差矩阵 Q
    std::deque<Matrix> QList(mesh->VertexNum(), Matrix(4, 4));
    for (int i = 0; i < mesh->VertexNum(); ++i)
    {
        vertex3f *v = mesh->m_vList[i];
        QList[i] = QEMatrix(v);

        // 对所有顶点标号，方便与 Q 对应
        v->identity = i;
    }

    // 获得网格上的所有边

    // 线性求解器
    LinearSolver solver;

    // 循环 times 次进行简化
    while (times > 0)
    {
        // 保存最优化的结果、误差和索引
        Matrix minQ(4, 4);
        Vector minv(4);
        double mine = FLT_MAX;
        int minInd = 0;

        // 循环半边，每个半边就是一对点
        for (int i = 0; i < mesh->m_eList.size(); i++)
        {
            half_edge *e = mesh->m_eList[i];
            vertex3f *v1 = e->v;
            vertex3f *v2 = destination(e);

            // 求解极小化问题的解 v，通过顶点标记获得对应的矩阵
            Matrix Q = QList[v1->identity] + QList[v2->identity];

            // 获得线性系统的矩阵并求解
            Matrix Q1 = Q;
            Q1(3, 0) = Q1(3, 1) = Q1(3, 2) = 0;
            Q1(3, 3) = 1;
            Vector b = {0, 0, 0, 1};
            Vector v = solver.solve(Q1, b, true);

            // 计算此时的误差，如果更小，就保存相关信息
            double err = dot(v, Q * v);
            if (err < mine)
            {
                minQ = Q;
                minv = v;
                mine = err;
                minInd = i;
            }
        }

        // 到这里，我们得到了具有最小误差的顶点合并，新顶点直接替换掉对应半边的起点
        // 只需要移动该顶点位置即可
        half_edge *e = mesh->m_eList[minInd];
        vertex3f *v1 = e->v;
        vertex3f *v2 = destination(e);
        v1->pos.x = minv[0];
        v1->pos.y = minv[1];
        v1->pos.z = minv[2];

        // 修改对应的误差矩阵
        QList[v1->identity] = minQ;

        // 获得需要操作的边，e11 与 e21 要配对，e12 与 e22 要配对
        half_edge *e11 = e->prev->twin;
        half_edge *e12 = (e->twin == nullptr) ? nullptr : e->twin->next->twin;
        half_edge *e21 = e->next->twin;
        half_edge *e22 = (e->twin == nullptr) ? nullptr : e->twin->prev->twin;

        // 进行删除操作，对于边收缩而言，会删除 1 到 2 个面
        // 第一个面是从 v1 指向 v2 的半边对应的面，这个面一定存在，因为 e 存在；第二个是该半边的反向边对应的面，不一定存在
        // 注意先删除 twin，不然 e 先被删除就成为 nullptr 了
        if (e->twin != nullptr)
            mesh->DeleteFace(e->twin->incFace);
        mesh->DeleteFace(e->incFace);

        // 获得相邻面，面上的 v2 全部改成 v1
        std::deque<face *> fList = adjacentFace(v2);
        for (int i = 0; i < fList.size(); ++i)
        {
            half_edge *tmpe = fList[i]->e;
            for (int j = 0; j < 3; j++)
            {
                if (tmpe->v == v2)
                {
                    tmpe->v = v1;
                    break;
                }
                tmpe = tmpe->next;
            }
        }

        // 令 v2 孤立
        v2->e = nullptr;

        // 配对这些半边
        if (e11 != nullptr && e21 != nullptr)
        {
            e11->twin = e21;
            e21->twin = e11;
        }
        if (e12 != nullptr && e22 != nullptr)
        {
            e12->twin = e22;
            e22->twin = e12;
        }

        // 检查半边个数，如果为零，则整个网格都被删除，返回错误
        if (mesh->m_eList.size() == 0)
            return false;

        times--;
    }

    // 最后要移除所有离散点
    mesh->DeleteDiscrete();
    return true;
}